library(ggplot2)
library(tidyverse)
library(patchwork)
library(lomb)
library(moments)
library(knitr)
library(oce)
library(fracdiff)
library(forecast)
library(astsa)
library(mFilter)
col <- "#2c3e50"
col2 <- "#18bc9c"
Zbiór danych pochodzi ze strony https://kp.gfz.de/en/data i zawiera historyczne pomiary aktywności geomagnetycznej Ziemi, która jest bezpośrednio skorelowana z występowaniem zórz polarnych. Analiza obejmuje okres od 1 stycznia 2018r. do końca 2025 roku. Jest to krótszy okres niż pełny cykl aktywności słonecznej (ok. 11 lat), jednak jest to zabieg celowy - w tej analizie skupimy się na częstszych cyklach, traktując pełny cykl aktywności słonecznej jako trend. Zestaw danych zawiera obserwacje o częstotliwości 3 godzinnej i obejmuje m.in.: Indeks Kp - indeks aktywności geomagnetycznej i Indeks Ap, który jest przeliczeniem indeksu Kp na skalę liniową. W tej analizie zajmiemy się średnią dzienną wartością Ap.
data <- read.csv("kp.csv", sep = ";", header = F)
dates <- as.Date(paste(data[,1], data[,2],data[,3], sep="-"))
ap <- data$V24
start_year <- as.numeric(format(dates[1], "%Y"))
start_day <- as.numeric(format(dates[1], "%j"))
xts <- ts(ap, frequency = 365, start = c(start_year, start_day))
ggplot(data.frame(Time = as.numeric(time(xts)), Ap = as.numeric(xts)), aes(x = Time, y = Ap)) +
geom_line(color = col) +
labs(title = "Aktywność geomagnetyczna (Indeks Ap)",
y = "Wartość indeksu Ap", x = "Czas") +
theme_bw()
stat_opisowe <- data.frame(
"Średnia" = mean(xts),
"Odch. stand." = sd(xts),
"Wsp. skośności" = skewness(xts),
"Kurtoza" = moments::kurtosis(xts),
"Min" = min(xts),
"Q1" = quantile(xts)[2],
"Mediana" = median(xts),
"Q3" = quantile(xts)[4],
"Max" = max(xts),
row.names = "Wartość"
)
kable(t(round(stat_opisowe,4)))
| Wartość | |
|---|---|
| Średnia | 9.1766 |
| Odch..stand. | 11.1310 |
| Wsp..skośności | 7.9004 |
| Kurtoza | 130.5495 |
| Min | 0.0000 |
| Q1 | 4.0000 |
| Mediana | 6.0000 |
| Q3 | 11.0000 |
| Max | 271.0000 |
p1 <- ggplot(data.frame(Time = as.numeric(time(xts)), Ap = as.numeric(xts)), aes(x = Ap)) +
geom_histogram(bins=100,fill = col ,color="white") +
labs(title = "Rozkład aktywności geomagnetycznej (Indeks Ap)",
y = "Częstość", x = "") +
theme_bw()
p2 <- ggplot(data.frame(Time = as.numeric(time(xts)), Ap = as.numeric(xts)), aes(y = "", x = Ap)) +
geom_boxplot(fill = "lightgrey", outlier.color = col, outlier.shape = 1) +
labs(, y="", x="Indeks Ap") +
theme_bw()
p1 / p2
Szereg nie jest stacjonarny, jego rozkład jest mocno skośny i istnieją obserwacje odstające. Zastosujemy logarytmizację w celu stabilizacji wariancji oraz usuniemy trend związany z pełnym cyklem aktywności słonecznej (trwającym 11 lat).
xts_log <- log(xts + 1)
model <- loess(xts_log ~ time(xts_log), span = 0.5)
xts_detrend <- residuals(model)
ggplot(data.frame(xts_log = xts_log, time = time(xts_log), fitted = model$fitted), aes(x = time)) +
geom_line(aes(y = xts_log), col="lightgrey") +
geom_line(aes(y = fitted), col=col, linewidth = 1) +
labs(title = "Szereg zlogarytmowany z dopasowanym trendem LOESS", x = "Czas", y = "") +
theme_bw()
mean(xts_detrend)
## [1] 0.005586405
Średnia po usunięciu trendu nie jest równa zero (funkcja LOESS nie działa jak klasyczna regresja, gdzie reszty są równe zero), dlatego stosujemy jeszcze centrowanie danych.
xts <- xts_detrend - mean(xts_detrend)
xts <- ts(xts[1:(length(xts) - 2)], frequency = 365, start = c(start_year, start_day))
ggplot(data.frame(xts = xts, time = time(xts)), aes(x = time, y = xts)) +
geom_line(col="lightgrey") +
labs(title = "Szereg scentrowany", x = "Czas", y = "") +
theme_bw()
n <- 2920
m <- floor(n/2)
fourier_log <- fft(xts_log[1:n])
fourier_log <- fourier_log[2:(m+1)]
y_log <- (1/(2*pi*n))* (Mod(fourier_log))^2
f <- (1:m)/n
#dane scentrowane
fourier <- fft(xts[1:2920])
fourier <- fourier[2:(m+1)]
y <- (1/(2*pi*n))* (Mod(fourier))^2
ggplot(data.frame(f = f, y=y, y_log = y_log), aes(x = f)) +
geom_line(aes(y = y_log, color = "Dane oryginalne"), linewidth = 0.7) +
geom_line(aes(y = y, color = "Dane scentrowane"), linewidth = 0.3) +
scale_color_manual(values = c("Dane oryginalne" = col,
"Dane scentrowane" = col2)) +
labs(title = "Periodogram: Dane oryginalne i scentrowane", x = "Czas", y = "", color = "") +
coord_cartesian(ylim = c(0, 13)) +
theme_bw() +
theme(legend.position = c(0.85, 0.85), legend.background = element_rect(fill = "white", color = "black", linewidth = 0.2))
ggplot(data.frame(f = f, y=y, y_log = y_log), aes(x = f)) +
geom_line(aes(y = y_log, color = "Dane oryginalne"), linewidth = 1) +
geom_line(aes(y = y, color = "Dane scentrowane"), linewidth = 1) +
scale_color_manual(values = c("Dane oryginalne" = col,
"Dane scentrowane" = col2)) +
labs(title = "Periodogram: Dane oryginalne i scentrowane (Zbliżenie)", x = "Czas", y = "", color = "") +
coord_cartesian(xlim = c(0, 0.02), ylim = c(0, 13)) +
theme_bw() +
theme(legend.position = c(0.85, 0.85), legend.background = element_rect(fill = "white", color = "black", linewidth = 0.2))
W przypadku danych oryginalnych (kolor granatowy), widmo jest całkowicie zdominowane przez energię skupioną przy częstotliwości zerowej.
n <- 2920
n2 <- 2900
m <- floor(n/2)
m2 <- floor(n2/2)
fourier_p <- fft(xts_log[1:n])
fourier_p <- fourier_p[2:(m+1)]
fourier_np <- fft(xts_log[1:n2])
fourier_np <- fourier_np[2:(m2+1)]
y_p <- (1/(2*pi*n))* (Mod(fourier_p))^2
y_np <- (1/(2*pi*n2))* (Mod(fourier_np))^2
f2 <- (1:m2)/n2
df_p <- data.frame(freq = f, spec = y_p, type = "N=2920 (podzielne)")
df_np <- data.frame(freq = f2, spec = y_np, type = "N=2900 (niepodzielne)")
df_combined <- rbind(df_p, df_np)
ggplot(df_combined, aes(x = freq, y = spec, color = type)) +
geom_line(linewidth = 0.7, alpha = 0.8) +
scale_color_manual(values = c("N=2920 (podzielne)" = col,
"N=2900 (niepodzielne)" = col2)) +
labs(title = "Periodogram w zależności od liczby obserwacji",
x = "Częstotliwość f", y = "Moc", color = "Liczba obserwacji N") +
theme_bw() +
theme(legend.position = c(0.8, 0.8),
legend.background = element_rect(fill = "white", color = "black", linewidth = 0.2))
ggplot(df_combined, aes(x = freq, y = spec, color = type)) +
geom_line(linewidth = 1, alpha = 0.8) +
scale_color_manual(values = c("N=2920 (podzielne)" = col,
"N=2900 (niepodzielne)" = col2)) +
coord_cartesian(xlim = c(0.0025, 0.01), ylim = c(0, 2.5)) +
labs(title = "Periodogram: Wyciek widma (Zbliżenie)",
x = "Częstotliwość f", y = "Moc", color = "Liczba obserwacji N") +
theme_bw() +
theme(legend.position = c(0.8, 0.8),
legend.background = element_rect(fill = "white", color = "black", linewidth = 0.2))
Skupiając się na trendzie półrocznym widzimy, że pik przy częstotliwości 0,0055 (182,5 dnia) jest węższy i wyższy dla N podzielnego przez 365. W przypadku N niepodzielnego przez 365 dochodzi do zjawiska wycieku widma, co widać na wykresie - pik jest niższy, a podstawa szersza. Analizując inne trendy nie obserwujemy tak silnego zjawiska, ponieważ są one naturalnie rozmyte (przykładowo aktywność Słońca trwa między 25 a 30 dni)
d <- 1/f
kable(slice_max(data.frame(Częstotliwość = f, Wartość = y, Dni = d), order_by = y, n=10))
| Częstotliwość | Wartość | Dni | |
|---|---|---|---|
| 17 | 0.0054795 | 2.2110420 | 182.500000 |
| 115 | 0.0390411 | 1.4271172 | 25.614035 |
| 98 | 0.0332192 | 1.0243959 | 30.103093 |
| 112 | 0.0380137 | 0.9180363 | 26.306306 |
| 100 | 0.0339041 | 0.8616870 | 29.494950 |
| 99 | 0.0335616 | 0.7755316 | 29.795918 |
| 103 | 0.0349315 | 0.7347768 | 28.627451 |
| 312 | 0.1065068 | 0.7215207 | 9.389067 |
| 105 | 0.0356164 | 0.6979422 | 28.076923 |
| 214 | 0.0729452 | 0.6638256 | 13.708920 |
Od teraz korzystamy tylko z danych scentrowanych po usunięciu trendu.
m <- floor(n/2)
f <- (1:m)/n
fourier <- fft(xts)
fourier <- fourier[2:(m+1)]
y <- (1/(2*pi*n))* (Mod(fourier))^2
cl <- (2 * y) / qchisq(1 - 0.025, df=2)
cu <- (2 * y) / qchisq(0.025, df=2)
plot_sur <- ggplot(data.frame(f = f, y=y), aes(x = f, y = y)) +
geom_line(col = col, linewidth = 0.8) +
labs(
x = "Częstotliwość",
y = " ",
title = "Periodogram naiwny"
) +
theme_bw()
ggplot(data.frame(f = f, y=y, cl = cl, cu = cu), aes(x = f)) +
geom_ribbon(aes(ymin = cl, ymax = cu), fill = "lightgrey") +
geom_line(aes(y = y), col = col, linewidth = 0.5) +
scale_y_log10() +
labs(
x = "Częstotliwość",
y = "Wartość (skala logarytmowana)",
title = "Periodogram naiwny z 95% przedziałem ufności"
) +
theme_bw()
p_welch <- pwelch(xts, plot = F)
p_welch$f <- p_welch$freq / 365
plot_welch <- ggplot(data.frame(f = p_welch$f, y=p_welch$spec), aes(x = f, y = y)) +
geom_line(col = col, linewidth = 0.8) +
labs(
x = "Częstotliwość",
y = " ",
title = "Periodogram metodą Welcha"
) +
theme_bw()
plot_welch
p_dan1 <- spec.pgram(
xts,
kernel = kernel("daniell", m = 5),
taper = 0,
detrend = F,
demean = F,
fast = F,
plot = F
)
p_dan2 <- spec.pgram(
xts,
kernel = kernel("daniell", m = c(3, 5)),
taper = 0,
detrend = F,
demean = F,
fast = F,
plot = F
)
p_dan1$f <- p_dan1$freq / 365
p_dan2$f <- p_dan2$freq / 365
plot_dan_1 <- ggplot(data.frame(f = p_dan1$f, y=p_dan1$spec), aes(x = f, y = y)) +
geom_line(col = col, linewidth = 0.8) +
labs(
x = "Częstotliwość",
y = " ",
title = "Periodogram wygładzony oknem Daniella (okno jednopoziomowe)"
) +
theme_bw()
plot_dan_2 <- ggplot(data.frame(f = p_dan2$f, y=p_dan2$spec), aes(x = f, y = y)) +
geom_line(col = col, linewidth = 0.8) +
labs(
x = "Częstotliwość",
y = " ",
title = "Periodogram wygładzony oknem Daniella (okno wielopoziomowe)"
) +
theme_bw()
plot_dan_1 / plot_dan_2
K <- 5
L <- floor(2920 / K)
data_matrix <- matrix(xts[1:(L*K)], nrow = L, ncol = K)
p_bartlett_temp <- spec.pgram(data_matrix[,1], taper=0, detrend = F, demean = F, plot=F)
avg_spec <- p_bartlett_temp$spec
for(i in 2:K) {
avg_spec <- avg_spec + spec.pgram(data_matrix[,i], taper=0, detrend = F, demean = F, plot=F)$spec
}
avg_spec <- avg_spec / K
p_bartlett <- p_bartlett_temp
p_bartlett$spec <- avg_spec
plot_bartlett <- ggplot(data.frame(f = p_bartlett$freq, y=p_bartlett$spec), aes(x = f, y = y)) +
geom_line(col = col, linewidth = 0.8) +
labs(
x = "Częstotliwość",
y = " ",
title = "Periodogram wygładzony przez rozłączne podokresy"
) +
theme_bw()
plot_bartlett
empty <- ggplot() + theme_void()
(plot_sur + plot_welch) / (plot_dan_1 + plot_dan_2) / (plot_bartlett + empty)
Periodogram naiwny wykazuje największą wariancję. Dzięki temu jego piki są bardzo wąskie, jednak zwiększa to ilość szumu przez co jest mniej czytelny. Periodogram wyznaczony metodą Welcha radzi sobie dobrze z usuwaniem szumu tła, nie rozmywając sygnału na boki. W przypadku periodogramów wygładzonych metodą Daniella piki są zdecydowanie szersze. Szczególnie widać to w przypadku okna wielopoziomowego. Dla periodogramu wygładzonego przez rozłączne podokresy szerokość piku jest wąska, natomiast wariancja jest zdecydowanie większa niż w przypadku Welcha i Daniella.
n <- length(xts)
idx <- sample(1:n, size = round(0.3 * n))
xts_NA <- xts
xts_NA[idx] <- NA
xts_NAomit <- na.omit(as.numeric(xts_NA))
p_naiwny_NA <- spec.pgram(
xts_NAomit,
log = "no",
detrend = FALSE,
demean = FALSE,
plot = FALSE,
)
p1_NA <- ggplot(data_frame(f = p_naiwny_NA$freq, s = p_naiwny_NA$spec), aes(f, s)) +
geom_line(col = col, linewidth = 0.6) +
labs(title = "Periodogram naiwny dla danych niekompletnych",
x = "Częstotliwość", y = " ") +
theme_bw()
ls <- lsp(xts_NA, type = "frequency", alpha = 0.05, plot = F, from = 0, to = 0.5)
p2_NA <- ggplot(data_frame(f = ls$scanned, s = ls$power), aes(f, s)) +
geom_line(col = col, linewidth = 0.6) +
geom_hline(yintercept = ls$sig.level, color = "red", linetype = "dashed") +
labs(title = "Periodogram Lomba–Scargle’a",
x = "Częstotliwość", y = " ") +
theme_bw()
p1_NA / p2_NA
Po usunięciu 30% obserwacji periodogram Lomba-Scargle’a poradził sobie lepiej niż periodogram naiwny. Przede wszystkim mimo braków w danych udało mu się odwzorować dwa najważniejsze piki związane z okresem półrocznym i trwającym ok. 27 dni. W przypadku periodogramu naiwnego, trend półroczny został odwzorowany poprawnie, jednak 27-dniowy jest rozmyty i zdecydowanie niższy niż powinien być.
t <- 1:300
x <- rnorm(300, 0, 0.1)
y <- sin(t/20) + x
v20 <- rep(1/20, 20)
v21 <- rep(1/21, 21)
plot(t, y, type = "l")
y_lag20 <- stats::filter(y, v20, sides=1)
lines(t, y_lag20, col="red")
y_lag21_centr <- stats::filter(y, v21, sides=2)
lines(t, y_lag21_centr, col="blue")
Średnie ruchome w skuteczny sposób usunęły szum i w odwzorowały trend okresowy. Średnia jednostrona składająca się z 20 obserwacji (kolor czerwony) wykazuje opóźnienie (przesunięcie fazowe) względem oryginału. W przypadku scentrowanej średniej ruchomej (kolor niebieski), jej wartości idealnie pokrywają się z fazą sygnału.
pkb <- read.csv("pkb.csv", sep=";")
plot(pkb$PKB, type = "l")
mv2 <- c(0.25, 0.5, 0.25)
y2 <- stats::filter(pkb$PKB, mv2, sides=2)
lines(y2, col="red")
mv3 <- c(1/3, 1/3, 1/3)
y3 <- stats::filter(pkb$PKB, mv3, sides=2)
lines(y3, col="blue")
mv4 <- c(1/8, 1/4, 1/4, 1/4, 1/8)
y4 <- stats::filter(pkb$PKB, mv4, sides=2)
lines(y4, col="green4")
mv5 <- c(1/5, 1/5, 1/5, 1/5, 1/5)
y5 <- stats::filter(pkb$PKB, mv5, sides=2)
lines(y5, col="orange")
Średnie ruchome rzędu 2 i 3 (czerwona i niebieska) wciąż podążają za sezonowymi wahaniami, nie eliminując ich całkowicie. Dopiero zastosowanie wag rzędu 4 (kolor zielony) pozwala na skuteczne odsezonowanie danych, co sugeruje cykl okresowości kwartalny. Jednocześnie przez wyższy rząd średniej ruchomej tracimy informację o dynamice zjawiska na końcach szeregu przez skracanie się wektora wyników. Średnia rzędu 5 nie radzi sobie najlepiej - pojawiają się sztuczne oscylacje w miejscu, w którym nie powinny.
ozon <- read.csv("ozon.csv", sep=" ")
x <- 1:200
y <- ozon[1:200,]
plot(x,y)
model <- lm(y ~ sin(2*x)+cos(2*x)+
sin(3*x)+cos(3*x)+
sin(4*x)+cos(4*x)+
sin(5*x)+cos(5*x)+
sin(6*x)+cos(6*x)+
sin(7*x)+cos(7*x)+
sin(8*x)+cos(8*x)+
sin(9*x)+cos(9*x))
y2 <- predict(model)
plot(x, y, type="l")
lines(y2, col="blue")
xnam <- paste0("sin(x*", 1:50,") + cos(x*", 1:50, ")")
fmla <- as.formula(paste("y ~ ", paste(xnam, collapse= "+")))
model2 <- lm(fmla)
y3 <- predict(model2)
lines(y3, col="red")
W przypadku dopasowania kombinacji harmonik sinusoidalnych o częstotliwościach k = 2,…9 (kolor niebieski) model dobrze wyłapuje ogólny rytm sezonowy jednak niedoszacowuje ekstremów. W przypadku modelu składającego się z 50 harmonik, model dopasowuje się zdecydowanie lepiej do danych, jednak niesie ze sobą ryzyko nadmiernego dopasowania (overfitting), przez co jakość prognozy poza próbą może być średniej jakości.
pkb <- read.csv("pkb.csv", sep=";")
n <- 45
pkb <- pkb[1:n,]
m <- floor(n/2)
fourier <- fft(pkb$PKB_detrended)
fourier <- fourier[2:(m+1)]
#periodogram
y <- (1/(2*pi*n))* (Mod(fourier))^2
f <- (1:m)/n
n2 <- 44
pkb2 <- pkb[1:n2,]
m2 <- floor(n2/2)
fourier2 <- fft(pkb2$PKB_detrended)
fourier2 <- fourier2[2:(m2+1)]
#periodogram
y2 <- (1/(2*pi*n2))* (Mod(fourier2))^2
f2 <- (1:m2)/n2
plot(f, y, type="l", col="blue", ylim= c(0, max(y2)))
lines(f2, y2, type="l", col="red")
abline(v=0.25, col="red", lty=3)
legend("topleft", legend=c("N=45 (niepodzielne)", "N=44 (podzielne)"),
col=c("blue", "red"), lty=1)
Dla N = 44, czyli podzielnego przez 4, pik sezonowy przypada dokładnie na częstotliwość 0,25. Dodatkowo pojawia się pik dla częstotliwości 0,5, który jest cyklem półrocznym, będącym wielokrotnością cyklu podstawowego. Badając dla N = 45, niepodzielnego przez 4, dochodzi do zjawisku wycieku widma. Widmo rozlewa się do sąsiednich częstotliwości przez co główny pik jest szerszy i niższy.
f[which.max(y)]
## [1] 0.2444444
f2[which.max(y2)]
## [1] 0.5
f2[which.max(y2[1:21])]
## [1] 0.25
sp500 <- read.csv("sp500.csv", sep=";")
sp500$RETURN <- as.numeric(sub("%", "", sp500$RETURN))
n <- length(sp500$RETURN)
m <- floor(n/2)
fourier <- fft(sp500$RETURN)
fourier <- fourier[2:(m+1)]
#periodogram
y <- (1/(2*pi*n))* (Mod(fourier))^2
f <- (1:m)/n
px <- ((1:6)/12)
py <- y[which(px==f)]
plot(f, y, type="l")
lines(px, py, type="p", col="red")
abline(v=px, col="red", lty=3)
Na periodogramie stóp zwrotu indeksu S&P 500 nie widać wyraźnych i dominujących pików. Częstotliwości 1/12, 2/12 itd. zaznaczone na czerwono nie wystają ponad przeciętny poziom szumu, co oznacza, że w danych miesięcznym składnik sezonowy jest nieistotny.
t <- seq(0, 0.296, by=0.001)
yt <- cos(2*pi*200*t) + rnorm(length(t), 0, 1)
yts <- ts(yt, frequency=1000)
yts <- yts - mean(yts)
#periodogram naiwny
p_naiwny <- spec.pgram(yts, detrend = T, taper= 0, plot=F)
#periodogram wygładzony oknem Daniella
p_daniell <- spec.pgram(yts, kernel("daniell", 3), taper = 0, detrend = F, plot = F)
#periodogram Welcha
p_welch <- pwelch(yts, noverlap = 10, plot = F)
plot(p_welch$freq, p_welch$spec, type ="l", ylim=c(0,0.02), col="red")
lines(p_naiwny$freq, p_naiwny$spec, col = "grey")
lines(p_daniell$freq, p_daniell$spec, col = "blue")
legend("topright", legend=c("Periodogram metodą Welcha", "Periodogram wygładzony oknem Daniella"),
col=c("red","blue"), lty=1)
Periodogram naiwny (kolor szary) charakteryzuje się największą wariancją. Periodogram wyznaczony metodą Welcha skutecznie wygładza szum tła i dodatkowo dobrze radzi sobie z pikiem. Periodogram wygładzony oknem Daniella charakteryzuje się szerokim i płaskim wierzchołkiem oraz bardzo stromymi krawędziami.
x <- fracdiff.sim(1000, ar = 0.1, ma = -0.4, d=0.4)$series
fit <- arfima(x)
tsdisplay(residuals(fit))
Reszty oscylują wokół zera. Brak istotnej autokorelacji.
Prognoza punktowa i przedziałowa
f_arfima <- forecast(fit, h=100)$mean
Arima
fit_arima <- auto.arima(x)
f_arima <- forecast(fit_arima, h=100)$mean
Porównanie
par(mfrow = c(1,2))
plot(forecast(fit, h=100))
plot(forecast(fit_arima, h=100))
par(mfrow = c(1,1))
Model ARFIMA z parametrem d = 0,4 wykazuje cechy długiej pamięci, a wpływ dawnych obserwacji wygasa powoli. Linia prognozy staje się płaska i widać na niej pamięć ostatnich obserwacji. Przedziały ufności są stabilne i rozszerzają się bardzo powoli. Natomiast model ARIMA szybko wraca do średniej, a przedziały ufności rozszerzają się gwałtownie.
df <- read.csv("infl_bezr.csv", sep=";")
xts <- ts(data = df$BEZR, frequency = 12)
xts <- xts - mean(xts)
cf <- mFilter::cffilter(xts, pl=3, pu=72)$cycle
cf2 <- mFilter::cffilter(xts, pl=6, pu=72)$cycle
hp <- mFilter::hpfilter(xts, type = "frequency", freq = 12)$cycle
p_raw <- spec.pgram(xts, taper = 0, pad=0, fast=F, detrend = F, plot=F)
p_cf <- spec.pgram(cf, taper = 0, pad=0, fast=F, detrend = F, plot=F)
p_cf2 <- spec.pgram(cf2, taper = 0, pad=0, fast=F, detrend = F, plot=F)
p_hp <- spec.pgram(hp, taper = 0, pad=0, fast=F, detrend = F, plot=F)
par(mfrow = c(2,2))
plot(p_raw, main="Periodogram: Szereg oryginalny")
plot(p_cf, main="Periodogram: Filtr CF (3-72)")
plot(p_cf2, main="Periodogram: Filtr CF (6-72)")
plot(p_hp, main="Periodogram: Filtr HP")
par(mfrow = c(1,1))
transfer_f_cf <- p_cf$spec / p_raw$spec
transfer_f_cf2 <- p_cf2$spec / p_raw$spec
transfer_f_hp <- p_hp$spec / p_raw$spec
freq <- p_raw$freq
plot(freq, transfer_f_cf, col="blue", type="l", lwd=2,
xlab="Częstotliwość", ylab=" ")
lines(freq, transfer_f_cf2, col="red", lwd=2)
lines(freq, transfer_f_hp, col="green", lwd=2)
legend("topleft", legend=c("CF (3-72)", "CF (6-72)", "HP"),
col=c("blue", "red", "green"), lty=1, lwd=2, bty="n")
Filtr CF (3-72) (kolor niebieski) przyjmuje najwyższe wartości dla częstotliwości do około 4, a następnie gwałtownie spada. Oznacza to, że filtr ten przepuszcza głównie średnie częstotliwości, a silnie tłumi zarówno trend jak i szum. W efekcie zachowuje jedynie wahania odpowiadające cyklom koniunkturalnym. Filtr CF (6-72) (kolor czerwony) przyjmuje wysokie wartości do częstotliwości około 2, po czym szybko zbliża się do zera. Oznacza, to że usuwa on jeszcze więcej krótkookresowowych wahań i koncentruje się tylko na wolniejszych cyklach. W przypadku filtra HP (kolor zielony) wartości rosną dla wyższych częstotliwości - filtr HP nie eliminuje szybkich wahań, a nawet częściowo je wzmacnia.
AirPassengers <- log(AirPassengers)
sarima(AirPassengers, 0, 1, 1, 0, 1, 1, 12)
## initial value -3.086228
## iter 2 value -3.267980
## iter 3 value -3.279950
## iter 4 value -3.285996
## iter 5 value -3.289332
## iter 6 value -3.289665
## iter 7 value -3.289672
## iter 8 value -3.289676
## iter 8 value -3.289676
## iter 8 value -3.289676
## final value -3.289676
## converged
## initial value -3.286464
## iter 2 value -3.286855
## iter 3 value -3.286872
## iter 4 value -3.286874
## iter 4 value -3.286874
## iter 4 value -3.286874
## final value -3.286874
## converged
## <><><><><><><><><><><><><><>
##
## Coefficients:
## Estimate SE t.value p.value
## ma1 -0.4018 0.0896 -4.4825 0
## sma1 -0.5569 0.0731 -7.6190 0
##
## sigma^2 estimated as 0.001348035 on 129 degrees of freedom
##
## AIC = -3.690069 AICc = -3.689354 BIC = -3.624225
##
sarima.for(xdata = AirPassengers, 12, 0, 1, 1, 0, 1, 1, 12)
## $pred
## Jan Feb Mar Apr May Jun Jul Aug
## 1961 6.110186 6.053775 6.171715 6.199300 6.232556 6.368779 6.507294 6.502906
## Sep Oct Nov Dec
## 1961 6.324698 6.209008 6.063487 6.168025
##
## $se
## Jan Feb Mar Apr May Jun
## 1961 0.03671562 0.04278291 0.04809072 0.05286830 0.05724856 0.06131670
## Jul Aug Sep Oct Nov Dec
## 1961 0.06513124 0.06873441 0.07215787 0.07542612 0.07855851 0.08157070
x <- diff(log(AirPassengers))
cf <- cffilter(x, pl = 3, pu = 72)$cycle
hp <- hpfilter(x, type = "frequency", freq = 12)$cycle
p_raw <- spec.pgram(x, taper = 0, pad = 0, fast = FALSE, detrend = FALSE, plot = FALSE)
p_cf <- spec.pgram(cf, taper = 0, pad = 0, fast = FALSE, detrend = FALSE, plot = FALSE)
p_hp <- spec.pgram(hp, taper = 0, pad = 0, fast = FALSE, detrend = FALSE, plot = FALSE)
plot(p_raw$freq, p_raw$spec, type="l", main="Widma: oryginał vs filtracja", lwd = 2)
lines(p_cf$freq, p_cf$spec, col="blue", lwd = 1)
lines(p_hp$freq, p_hp$spec, col="green", lwd = 1)
legend("topright",
legend=c("Oryginał", "CF", "HP"),
col=c("black","blue","green"),
lwd = 1, bty="n")
Podobnie jak w przypadku poprzedniego zadania filtr CF zachowuje wybrane składowe cykliczne, natomiast HP przepuszcza znaczną częśc szumu o wysokiej częstotliwości.
transfer_cf <- p_cf$spec / p_raw$spec
transfer_hp <- p_hp$spec / p_raw$spec
freq <- p_raw$freq
plot(freq, transfer_hp, type="l", col="blue", lwd=2,
xlab="Częstotliwość", ylab="Funkcja przenoszenia mocy")
lines(freq, transfer_cf, col="green", lwd=2)
legend("topright",
legend=c("HP", "CF"),
col=c("blue","green"),
lty=1, lwd=2, bty="n")
ker <- kernel("modified.daniell", 6)
p_smooth <- kernapply(p_raw$spec, ker)
freq_s <- freq[1:length(p_smooth)]
plot(freq_s, p_smooth, type="l",
main="Wygładzony periodogram (Daniell)",
xlab="Częstotliwość", ylab="Widmo")
Dla niskich częstotliwości wartości widma są stosunkowo wysokie. Oznacza to, że w danych występują istotne wolne wahania o długim okresie. Następnie widoczny jest wyraźny spadek wartości widma, co wskazuje na silne tłumienie średnich częstotliwości. Dla wyższych częstotliwości widmo ponownie rośnie, jednak jego poziom pozostaje znacznie niższy niż osiągnięty dla niskich częstotliwości.
dane <- read.csv("infl_bezr.csv", sep=";")
infl_ts <- ts(dane$INFL, frequency = 12, start = c(1, 1))
bezr_ts <- ts(dane$BEZR, frequency = 12, start = c(1, 1))
infl_d <- diff(infl_ts)
bezr_d <- diff(bezr_ts)
xts <- ts.intersect(infl_d, bezr_d)
t <- time(infl_d)
p_raw <- spec.pgram(infl_d, plot = FALSE)
p6 <- spec.pgram(infl_d, spans = 6, plot = FALSE)
p8 <- spec.pgram(infl_d, spans = 8, plot = FALSE)
omega <- 2 * pi * p_raw$freq
plot(omega, p_raw$spec, type="l", col="gray", lwd=1,
main="Periodogram inflacji",
xlab="Częstotliwość (radiany)",
ylab="Gęstość widmowa")
lines(2*pi*p6$freq, p6$spec, col="blue", lwd=2)
lines(2*pi*p8$freq, p8$spec, col="red", lwd=2)
legend("topright",
legend=c("surowy", "span=6", "span=8"),
col=c("gray","blue","red"),
lwd=c(1,2,2))
p <- spec.pgram(xts, spans = 6, taper=0, pad=0, fast=FALSE, plot = F)
omega <- 2 * pi * p$freq
plot(omega, p$coh, type="l",
xlab="Częstotliwość (radiany)",
ylab="Koherencja",
main="Koherencja inflacja–bezrobocie")
Koherencja przyjmuje najwyższe wartości dla niskich częstotliwości. Oznacza to, że inflacja i bezrobocie są umiarkowanie powiązane w długim okresie.
idx <- which.max(p$coh)
phi <- p$phase[idx]
omega <- 2 * pi * p$f[idx]
tau <- phi / omega
tau
## [1] 0.6013668
Opóźnienie czasowe (tau) wynosi 0,6 roku, czyli około 7 miesięcy - bezrobocie reaguje z opóźnieniem na zmiany inflacji.